#!/usr/bin/python
# coding: utf-8

import sys
from math import sqrt
from optparse import Option, OptionParser

from sflib.py_finder import create_norm_mx, create_mx, print_mx, inf_mx
from sflib.cpp_finder import DPResultExt
from sflib.mx_detect import iter_base_seqs, uniq
from shuffler import pshuffle

ALPHABET = 'actg'
PERIOD   = 3

def extend_backtrace(backtrace, win_length):
    last_shift = backtrace[-1]
    backtrace.extend([last_shift] * ( ((len(backtrace) + win_length)/PERIOD)*PERIOD - len(backtrace) ) )


def find_random(options, (seq_id, seq)):
    max_reward = 0
    max_Z = 0
    max_I = 0
    best_base_seq = None
    best_res = DPResultExt()
    for base_seq in iter_base_seqs(ALPHABET, PERIOD, options.win_length, options.mx_parts_num):
        #print base_seq
        #print seq
        #print 

        #from sflib.py_finder import create_dp_mx, read_dp_mx
        #dp_mx = create_dp_mx(options, base_seq, seq)
        ##print_dp_mx(dp_mx)
        #reward, start_pos, shifts = read_dp_mx(dp_mx)
        #cnt_shifts = max(0, len(shifts) - 1)
        #subseq_len = sum([i[1] for i in shifts])
        #if options.verbose:
        #    print "CURRENT: %3i\t%7.1f\t%s\t%i\t%i" % (cnt_shifts, reward, "[%s:%s]" % (start_pos, start_pos+subseq_len), subseq_len, shifts[0][1])
        ##subseq = seq[start_pos:start_pos+subseq_len]
        ##print subseq_len, backtrace
        ##print subseq
        ##print_res(subseq, shifts)

        from sflib.cpp_finder import create_dp_mx
        res = create_dp_mx(options, base_seq, seq, local=True)
        reward = res.reward
        cnt_shifts = res.sh_num
        subseq = seq[res.left_x:res.right_x+options.win_length+1]
        #print len(subseq)
        mk_res = create_dp_mx(options, base_seq, subseq, local=False)
        #print mk_res.reward, mk_res.left_x, mk_res.right_x
        # Monte-Karlo
        MK_NUM = 50
        rewards = []
        for i in range(MK_NUM):
            rseq = pshuffle(PERIOD, subseq)
            mk_res = create_dp_mx(options, base_seq, rseq, local=False)
            rewards.append( mk_res.reward )
        M = sum(rewards) / MK_NUM
        D = sum((reward - M)**2 for reward in rewards ) / MK_NUM
        Z = (res.reward - M) / sqrt(D)
        sub_mx = create_mx(subseq)
        I = 2*sqrt(inf_mx(sub_mx)) - sqrt(2*6-1)

        if options.verbose:
            print "CURRENT: %3i\t%7.1f\t%s\t%i\t%i\t%3.2f\t%3.2f" % (cnt_shifts, reward, "[%s:%s]" % (res.left_x, res.right_x), res.right_x-res.left_x, res.sh1, Z, I)
        
        if reward > max_reward:
            max_Z, max_I = Z, I
            max_reward = reward
            best_base_seq = base_seq
            best_res = res
    print "\t".join( map(str, (seq_id, len(seq), best_res.sh_num, best_res.sh1, best_res.left_x, best_res.right_x, best_res.reward, max_Z, max_I) ))
    if best_base_seq:
        print best_base_seq
        print_mx(create_norm_mx(best_base_seq))


if __name__ == '__main__':
    parser = OptionParser(usage='%prog [-r REWARD -p PENALTY -l WIN_LENGTH -s SEQUENCE]')
    parser.add_option(Option(
        '-r', '--reward',
        dest = 'reward',
        type = int,
        default = 1,
    ))
    parser.add_option(Option(
        '-p', '--penalty',
        dest = 'penalty',
        type = int,
        default = 3,
    ))
    parser.add_option(Option(
        '-l', '--window-length',
        dest = 'win_length',
        type = int,
        default = 36,
    ))
    parser.add_option(Option(
        '-s', '--sequence',
        dest = 'sequence',
    ))
    parser.add_option(Option(
        '-m', '--mx-parts',
        dest = 'mx_parts_num',
        type = int,
    ))
    parser.add_option(Option(
        '-v', '--verbose',
        dest = 'verbose',
        action = 'store_true',
        default = False,
    ))

    opts, args = parser.parse_args()
    if opts.win_length % PERIOD != 0:
        print >> sys.stderr, "WIN_LENGTH is not divided by %s" % PERIOD
        sys.exit(1)
        
    if opts.sequence:
        find_random(opts, opts.sequence)
    else:
        for line in sys.stdin:
            id, line = line.strip().split()
            find_random(opts, (id, line))
